HerbPt1 16S Figures

```r
require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)
require(dplyr)
require(ggplot2)
require(microbiome)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBtaWNyb2Jpb21lXG5XYXJuaW5nIGluIGxpYnJhcnkocGFja2FnZSwgbGliLmxvYyA9IGxpYi5sb2MsIGNoYXJhY3Rlci5vbmx5ID0gVFJVRSwgbG9naWNhbC5yZXR1cm4gPSBUUlVFLCAgOlxuICB0aGVyZSBpcyBubyBwYWNrYWdlIGNhbGxlZCDigJhtaWNyb2Jpb21l4oCZXG4ifQ== -->

Loading required package: microbiome Warning in library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, : there is no package called ‘microbiome’

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucmVxdWlyZSh2ZWdhbilcbmBgYFxuYGBgIn0= -->

```r
```r
require(vegan)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Load data

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucHNfZG1uIDwtIHJlYWRSRFMoXCIvVXNlcnMvZ29yZG9uY3VzdGVyL0Rlc2t0b3AvR2l0X1Byb2plY3RzL0hlcmJpY2lkZV9NaWNyb2Jlc19QVDEvZGF0YS9QaHlsb3NlcU9iamVjdHMvSVRTL0RNTl9lc3RzX0lUUy5SZGF0YVwiKVxuYGBgIn0= -->

```r
ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/DMN_ests_ITS.Rdata")

Warning message:
In rm(list = objects, envir = env) : object '*tmp*' not found

sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_dmn)$herb_time<-paste(sample_data(ps_dmn)$Herbicide, sample_data(ps_dmn)$Time, sep = "_")

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_dmn)$Mode<-sample_data(ps_dmn)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_dmn)$Mode<- as.factor(values[match(sample_data(ps_dmn)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_dmn)$Herbicide <- as.factor(values[match(sample_data(ps_dmn)$Herbicide, index)])


ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_rare_ITS.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_rare)$herb_time<-paste(sample_data(ps_rare)$Herbicide, sample_data(ps_rare)$Time, sep = "_")


#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Herbicide <- as.factor(values[match(sample_data(ps_rare)$Herbicide, index)])


ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/ITS/HerbPt1_hel_trans_ITS.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Aatrex", "Clarity", "Hand","Non-Treated","Roundup Powermax"))
sample_data(ps_trans)$herb_time<-paste(sample_data(ps_trans)$Herbicide, sample_data(ps_trans)$Time, sep = "_")


#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_trans)$Mode<-sample_data(ps_trans)$Herbicide

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Hand", "Non-Treated")

sample_data(ps_trans)$Mode<- as.factor(values[match(sample_data(ps_trans)$Mode, index)])

index <- c("Clarity", "Roundup Powermax", "Aatrex", "Hand", "Non-Treated")
values <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")

sample_data(ps_trans)$Herbicide <- as.factor(values[match(sample_data(ps_trans)$Herbicide, index)])

Remove samples that are outlines or under sequenced.

ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G009SG")
ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G095SG")
ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G123SG")
ps_dmn <-  subset_samples(ps_dmn, sample_names(ps_dmn) != "G129SG")

ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G009SG")
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G095SG")
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G123SG")
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G129SG")
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)

ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transformed to successfully ordinate. At this point, the transformed data will not ordinate. This section is full dataset ordinations.

ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)

Run 0 stress 0.1637061 
Run 1 stress 0.1638069 
... Procrustes: rmse 0.04511874  max resid 0.2959219 
Run 2 stress 0.1639952 
... Procrustes: rmse 0.04572109  max resid 0.2996605 
Run 3 stress 0.1636993 
... New best solution
... Procrustes: rmse 0.003124  max resid 0.03552253 
Run 4 stress 0.1636321 
... New best solution
... Procrustes: rmse 0.01768586  max resid 0.1058303 
Run 5 stress 0.1634893 
... New best solution
... Procrustes: rmse 0.01749345  max resid 0.1018927 
Run 6 stress 0.1636166 
... Procrustes: rmse 0.04806717  max resid 0.3003009 
Run 7 stress 0.1633705 
... New best solution
... Procrustes: rmse 0.04470724  max resid 0.2996104 
Run 8 stress 0.1631225 
... New best solution
... Procrustes: rmse 0.01730151  max resid 0.1523431 
Run 9 stress 0.1636991 
Run 10 stress 0.1640706 
Run 11 stress 0.1635635 
... Procrustes: rmse 0.01665523  max resid 0.203607 
Run 12 stress 0.163487 
... Procrustes: rmse 0.04938977  max resid 0.2995861 
Run 13 stress 0.1642535 
Run 14 stress 0.1633718 
... Procrustes: rmse 0.01725707  max resid 0.1522712 
Run 15 stress 0.1640259 
Run 16 stress 0.1646096 
Run 17 stress 0.1634884 
... Procrustes: rmse 0.04947432  max resid 0.2997511 
Run 18 stress 0.1631224 
... New best solution
... Procrustes: rmse 0.001164614  max resid 0.007901557 
... Similar to previous best
Run 19 stress 0.1631223 
... New best solution
... Procrustes: rmse 0.000894574  max resid 0.007191738 
... Similar to previous best
Run 20 stress 0.1632532 
... Procrustes: rmse 0.007653146  max resid 0.08811271 
*** Solution reached

ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1837394 
Run 1 stress 0.1837536 
... Procrustes: rmse 0.004222286  max resid 0.03854408 
Run 2 stress 0.1859057 
Run 3 stress 0.1837379 
... New best solution
... Procrustes: rmse 0.0006553106  max resid 0.004505322 
... Similar to previous best
Run 4 stress 0.1856626 
Run 5 stress 0.1861936 
Run 6 stress 0.1862078 
Run 7 stress 0.1859133 
Run 8 stress 0.1837401 
... Procrustes: rmse 0.0007329372  max resid 0.004394223 
... Similar to previous best
Run 9 stress 0.1837556 
... Procrustes: rmse 0.004232884  max resid 0.03808415 
Run 10 stress 0.1837291 
... New best solution
... Procrustes: rmse 0.004668855  max resid 0.03839136 
Run 11 stress 0.1862108 
Run 12 stress 0.1856156 
Run 13 stress 0.1837321 
... Procrustes: rmse 0.001869121  max resid 0.0172899 
Run 14 stress 0.185933 
Run 15 stress 0.1859134 
Run 16 stress 0.1837366 
... Procrustes: rmse 0.003934702  max resid 0.03800597 
Run 17 stress 0.183736 
... Procrustes: rmse 0.004072096  max resid 0.03782902 
Run 18 stress 0.1862073 
Run 19 stress 0.1837283 
... New best solution
... Procrustes: rmse 0.003314344  max resid 0.02283046 
Run 20 stress 0.1861278 
Run 21 stress 0.1837368 
... Procrustes: rmse 0.004202959  max resid 0.03676119 
Run 22 stress 0.1862196 
Run 23 stress 0.1837526 
... Procrustes: rmse 0.002014267  max resid 0.0209574 
Run 24 stress 0.1837541 
... Procrustes: rmse 0.001993472  max resid 0.02133405 
Run 25 stress 0.183726 
... New best solution
... Procrustes: rmse 0.0003499443  max resid 0.002495313 
... Similar to previous best
*** Solution reached

ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.01, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", k=3, trymax= 300, maxit = 1000)

Run 0 stress 0.1050141 
Run 1 stress 0.103327 
... New best solution
... Procrustes: rmse 0.03564627  max resid 0.1241836 
Run 2 stress 0.1033284 
... Procrustes: rmse 0.0004026237  max resid 0.003428707 
... Similar to previous best
Run 3 stress 0.1033322 
... Procrustes: rmse 0.000709747  max resid 0.004562807 
... Similar to previous best
Run 4 stress 0.1034339 
... Procrustes: rmse 0.00953169  max resid 0.05735493 
Run 5 stress 0.103324 
... New best solution
... Procrustes: rmse 0.002440901  max resid 0.0170954 
Run 6 stress 0.1034354 
... Procrustes: rmse 0.007967689  max resid 0.05434794 
Run 7 stress 0.1033282 
... Procrustes: rmse 0.00265007  max resid 0.01838042 
Run 8 stress 0.1034361 
... Procrustes: rmse 0.007981759  max resid 0.05521913 
Run 9 stress 0.1034359 
... Procrustes: rmse 0.008040017  max resid 0.05456366 
Run 10 stress 0.1033249 
... Procrustes: rmse 0.002039597  max resid 0.01573552 
Run 11 stress 0.1050151 
Run 12 stress 0.1034358 
... Procrustes: rmse 0.00800128  max resid 0.05441324 
Run 13 stress 0.1050135 
Run 14 stress 0.1068755 
Run 15 stress 0.1050132 
Run 16 stress 0.1033314 
... Procrustes: rmse 0.003020749  max resid 0.02047583 
Run 17 stress 0.1050126 
Run 18 stress 0.103327 
... Procrustes: rmse 0.002759821  max resid 0.02293387 
Run 19 stress 0.1033299 
... Procrustes: rmse 0.002844598  max resid 0.01951695 
Run 20 stress 0.103436 
... Procrustes: rmse 0.008047048  max resid 0.05456818 
Run 21 stress 0.1050133 
Run 22 stress 0.1050145 
Run 23 stress 0.1034335 
... Procrustes: rmse 0.007696545  max resid 0.05320767 
Run 24 stress 0.1050135 
Run 25 stress 0.1034323 
... Procrustes: rmse 0.007414607  max resid 0.05157759 
Run 26 stress 0.1033298 
... Procrustes: rmse 0.002836647  max resid 0.01947178 
Run 27 stress 0.1050129 
Run 28 stress 0.1050133 
Run 29 stress 0.1050146 
Run 30 stress 0.1033263 
... Procrustes: rmse 0.002509226  max resid 0.02018023 
Run 31 stress 0.103322 
... New best solution
... Procrustes: rmse 0.00036451  max resid 0.002650074 
... Similar to previous best
*** Solution reached

create alphadiversity tables

alpha_div_md$Herbicide 

  [1] Dicamba             Dicamba             Dicamba             Handweeded          Handweeded         
  [6] Handweeded          Non-Treated         Non-Treated         Atrazine-Mesotrione Atrazine-Mesotrione
 [11] Glyphosate          Glyphosate          Glyphosate          Glyphosate          Glyphosate         
 [16] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded         
 [21] Dicamba             Dicamba             Dicamba             Non-Treated         Non-Treated        
 [26] Non-Treated         Non-Treated         Non-Treated         Non-Treated         Glyphosate         
 [31] Glyphosate          Glyphosate          Dicamba             Dicamba             Dicamba            
 [36] Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded          Handweeded         
 [41] Glyphosate          Glyphosate          Glyphosate          Handweeded          Handweeded         
 [46] Non-Treated         Non-Treated         Non-Treated         Dicamba             Dicamba            
 [51] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Dicamba             Dicamba            
 [56] Handweeded          Handweeded          Handweeded          Non-Treated         Non-Treated        
 [61] Non-Treated         Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate         
 [66] Glyphosate          Glyphosate          Glyphosate          Glyphosate          Glyphosate         
 [71] Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded         
 [76] Handweeded          Dicamba             Dicamba             Dicamba             Non-Treated        
 [81] Non-Treated         Non-Treated         Non-Treated         Glyphosate          Glyphosate         
 [86] Dicamba             Dicamba             Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione
 [91] Handweeded          Handweeded          Handweeded          Glyphosate          Glyphosate         
 [96] Glyphosate          Handweeded          Handweeded          Non-Treated         Non-Treated        
[101] Non-Treated         Dicamba             Dicamba             Atrazine-Mesotrione Dicamba            
[106] Dicamba             Handweeded          Handweeded          Handweeded          Non-Treated        
[111] Non-Treated         Atrazine-Mesotrione Atrazine-Mesotrione Atrazine-Mesotrione Glyphosate         
[116] Glyphosate          Glyphosate          Glyphosate          Glyphosate          Atrazine-Mesotrione
[121] Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded          Handweeded         
[126] Dicamba             Dicamba             Dicamba             Non-Treated         Non-Treated        
[131] Non-Treated         Non-Treated         Non-Treated         Glyphosate          Glyphosate         
[136] Glyphosate          Dicamba             Dicamba             Dicamba             Atrazine-Mesotrione
[141] Atrazine-Mesotrione Atrazine-Mesotrione Handweeded          Handweeded          Glyphosate         
[146] Glyphosate          Glyphosate          Handweeded          Handweeded          Handweeded         
[151] Non-Treated         Non-Treated         Non-Treated         Dicamba             Atrazine-Mesotrione
[156] Atrazine-Mesotrione Atrazine-Mesotrione
Levels: Non-Treated Handweeded Atrazine-Mesotrione Dicamba Glyphosate

Shannon Div plots - no significant differences among herbicide treatments at any of the three time points

ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Shannon.pdf")

Saving 7.29 x 4.51 in image

aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)

summary(aov_t1)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4  0.389 0.09720    1.26  0.299
Residuals   48  3.704 0.07716               

aov_t2<-aov(Shannon~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)

summary(aov_t2)

            Df Sum Sq Mean Sq F value Pr(>F)
Herbicide    4 0.1234 0.03086   0.508   0.73
Residuals   46 2.7936 0.06073               

aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)

summary(aov_t3)

            Df Sum Sq Mean Sq F value Pr(>F)  
Herbicide    4 0.5376 0.13440   2.551  0.051 .
Residuals   48 2.5287 0.05268                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(aov_t3, "Herbicide")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3", ])

$Herbicide
                                        diff         lwr       upr     p adj
Handweeded-Non-Treated           0.086844630 -0.19738244 0.3710717 0.9078666
Atrazine-Mesotrione-Non-Treated  0.129708505 -0.14882205 0.4082391 0.6803851
Dicamba-Non-Treated             -0.002861275 -0.30174872 0.2960262 0.9999999
Glyphosate-Non-Treated           0.273363070 -0.01086400 0.5575901 0.0647237
Atrazine-Mesotrione-Handweeded   0.042863876 -0.22867317 0.3144009 0.9914461
Dicamba-Handweeded              -0.089705904 -0.38208716 0.2026754 0.9066040
Glyphosate-Handweeded            0.186518440 -0.09085878 0.4638957 0.3282566
Dicamba-Atrazine-Mesotrione     -0.132569780 -0.41941651 0.1542769 0.6864852
Glyphosate-Atrazine-Mesotrione   0.143654565 -0.12788248 0.4151916 0.5679081
Glyphosate-Dicamba               0.276224345 -0.01615691 0.5686056 0.0724032

Adonis testing of herbicide treatments by time point

ps_adonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
  print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
  vegan::adonis(physeq_dist ~ Herbicide * Time + Total_Weed_Veg , data = md_tab, permutations = 1000)
}

#ps_adonis(ps_rare_sub)
#remove one sample with no vegetation measurement. 
ps_rare_sub_57<-subset_samples(ps_rare_sub, sample_names(ps_rare_sub) != "G065SG")
ps_adonis(ps_rare_sub_57)
#ps_adonis(ps_trans_sub)
ps_dmn_57<-subset_samples(ps_dmn, sample_names(ps_dmn) != "G065SG")
#ps_adonis(ps_dmn)
ps_adonis(ps_dmn_57)

Ordination plots DMN by time point

ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1335559 
Run 1 stress 0.1309245 
... New best solution
... Procrustes: rmse 0.06976384  max resid 0.2060325 
Run 2 stress 0.1315098 
Run 3 stress 0.132755 
Run 4 stress 0.134753 
Run 5 stress 0.1308758 
... New best solution
... Procrustes: rmse 0.0538814  max resid 0.2344595 
Run 6 stress 0.1331698 
Run 7 stress 0.1323904 
Run 8 stress 0.1346675 
Run 9 stress 0.1350022 
Run 10 stress 0.1322358 
Run 11 stress 0.1334733 
Run 12 stress 0.1315248 
Run 13 stress 0.1331407 
Run 14 stress 0.1327605 
Run 15 stress 0.1332143 
Run 16 stress 0.1369199 
Run 17 stress 0.1339829 
Run 18 stress 0.1331858 
Run 19 stress 0.1360529 
Run 20 stress 0.1352334 
Run 21 stress 0.1308759 
... Procrustes: rmse 0.001713885  max resid 0.006848176 
... Similar to previous best
*** Solution reached

T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_dmn$plot + theme_classic()

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T1.pdf")

Saving 7.29 x 4.51 in image

ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1325262 
Run 1 stress 0.1325269 
... Procrustes: rmse 0.0001074309  max resid 0.000408716 
... Similar to previous best
Run 2 stress 0.1325262 
... New best solution
... Procrustes: rmse 0.001015874  max resid 0.005488789 
... Similar to previous best
Run 3 stress 0.1327962 
... Procrustes: rmse 0.01038273  max resid 0.06195753 
Run 4 stress 0.1338926 
Run 5 stress 0.1392078 
Run 6 stress 0.1338929 
Run 7 stress 0.1325264 
... Procrustes: rmse 0.001066291  max resid 0.005401677 
... Similar to previous best
Run 8 stress 0.1325261 
... New best solution
... Procrustes: rmse 0.0001868934  max resid 0.001066974 
... Similar to previous best
Run 9 stress 0.1334238 
Run 10 stress 0.1325259 
... New best solution
... Procrustes: rmse 0.0007505591  max resid 0.004597761 
... Similar to previous best
Run 11 stress 0.132925 
... Procrustes: rmse 0.01264762  max resid 0.07728014 
Run 12 stress 0.1438152 
Run 13 stress 0.132526 
... Procrustes: rmse 0.0003451055  max resid 0.001743723 
... Similar to previous best
Run 14 stress 0.1325259 
... Procrustes: rmse 0.0003370863  max resid 0.00167214 
... Similar to previous best
Run 15 stress 0.1334866 
Run 16 stress 0.1325266 
... Procrustes: rmse 0.0009514394  max resid 0.005508805 
... Similar to previous best
Run 17 stress 0.1325265 
... Procrustes: rmse 0.0008955128  max resid 0.005161325 
... Similar to previous best
Run 18 stress 0.1449845 
Run 19 stress 0.1342762 
Run 20 stress 0.1325263 
... Procrustes: rmse 0.0003284059  max resid 0.001374251 
... Similar to previous best
*** Solution reached

T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T2.pdf")

Saving 7.29 x 4.51 in image

#this time point needs to be checked out. The ordination is not working properly. 
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Run 0 stress 0.1451733 
Run 1 stress 0.1459287 
Run 2 stress 0.1510555 
Run 3 stress 0.1501906 
Run 4 stress 0.1503806 
Run 5 stress 0.1459215 
Run 6 stress 0.1528692 
Run 7 stress 0.1451733 
... Procrustes: rmse 0.00281782  max resid 0.01088482 
Run 8 stress 0.1502075 
Run 9 stress 0.1507344 
Run 10 stress 0.1478567 
Run 11 stress 0.1490906 
Run 12 stress 0.1451724 
... New best solution
... Procrustes: rmse 0.0002730257  max resid 0.001361369 
... Similar to previous best
Run 13 stress 0.1507808 
Run 14 stress 0.146477 
Run 15 stress 0.1510089 
Run 16 stress 0.1507665 
Run 17 stress 0.1518463 
Run 18 stress 0.1459425 
Run 19 stress 0.1451714 
... New best solution
... Procrustes: rmse 0.001928599  max resid 0.006972218 
... Similar to previous best
Run 20 stress 0.1459215 
*** Solution reached

T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_dmn_T3.pdf")

Saving 7.29 x 4.51 in image

Ordination plots on rarefied data by time point.

ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.15407 
Run 1 stress 0.1547275 
Run 2 stress 0.1560278 
Run 3 stress 0.154304 
... Procrustes: rmse 0.01759432  max resid 0.08547907 
Run 4 stress 0.1552432 
Run 5 stress 0.1559291 
Run 6 stress 0.1560286 
Run 7 stress 0.1560271 
Run 8 stress 0.1540703 
... Procrustes: rmse 0.0005180512  max resid 0.002163285 
... Similar to previous best
Run 9 stress 0.1548105 
Run 10 stress 0.1540705 
... Procrustes: rmse 0.0005653613  max resid 0.002608133 
... Similar to previous best
Run 11 stress 0.1581519 
Run 12 stress 0.1559287 
Run 13 stress 0.1559305 
Run 14 stress 0.1595564 
Run 15 stress 0.1540698 
... New best solution
... Procrustes: rmse 0.0002352417  max resid 0.000845193 
... Similar to previous best
Run 16 stress 0.1585595 
Run 17 stress 0.1540701 
... Procrustes: rmse 0.0004692716  max resid 0.002323505 
... Similar to previous best
Run 18 stress 0.1541024 
... Procrustes: rmse 0.004124925  max resid 0.02278237 
Run 19 stress 0.1543293 
... Procrustes: rmse 0.008122859  max resid 0.02901827 
Run 20 stress 0.158152 
*** Solution reached

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare$plot + theme_classic()

T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T1_rare_plot<-T1_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5)  + guides(color=guide_legend("Treatment")) + xlab("NMDS 1") + ylab("NMDS 2")
T1_rare_plot

Warning: Removed 3 rows containing missing values (geom_point).

library(cowplot)
my_legend <- get_legend(T1_rare_plot)

Warning: Removed 3 rows containing missing values (geom_point).

library(ggpubr)
as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordinationlegend.pdf")

Saving 7.29 x 4.51 in image

T1_rare_plot<-T1_rare_plot + theme(legend.position = "none")
T1_rare_plot

Warning: Removed 3 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 3 rows containing missing values (geom_point).

ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare_sub, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1716481 
Run 1 stress 0.1712224 
... New best solution
... Procrustes: rmse 0.07764546  max resid 0.2246421 
Run 2 stress 0.1706262 
... New best solution
... Procrustes: rmse 0.07905606  max resid 0.2661682 
Run 3 stress 0.1713448 
Run 4 stress 0.1709004 
... Procrustes: rmse 0.05656489  max resid 0.2944112 
Run 5 stress 0.1702319 
... New best solution
... Procrustes: rmse 0.06383277  max resid 0.2810499 
Run 6 stress 0.1706227 
... Procrustes: rmse 0.06377346  max resid 0.2806419 
Run 7 stress 0.1715869 
Run 8 stress 0.1713431 
Run 9 stress 0.1760045 
Run 10 stress 0.1710628 
Run 11 stress 0.1713372 
Run 12 stress 0.1709006 
Run 13 stress 0.1730992 
Run 14 stress 0.1721699 
Run 15 stress 0.1720581 
Run 16 stress 0.1702893 
... Procrustes: rmse 0.05501828  max resid 0.2833975 
Run 17 stress 0.1724195 
Run 18 stress 0.1707193 
... Procrustes: rmse 0.01192387  max resid 0.06232137 
Run 19 stress 0.170143 
... New best solution
... Procrustes: rmse 0.006859545  max resid 0.02883077 
Run 20 stress 0.172169 
Run 21 stress 0.1703025 
... Procrustes: rmse 0.05780913  max resid 0.2825995 
Run 22 stress 0.1709835 
Run 23 stress 0.1702878 
... Procrustes: rmse 0.05725869  max resid 0.2827813 
Run 24 stress 0.1706262 
... Procrustes: rmse 0.06561008  max resid 0.2790836 
Run 25 stress 0.1710053 
Run 26 stress 0.1715864 
Run 27 stress 0.1770233 
Run 28 stress 0.1728199 
Run 29 stress 0.1702891 
... Procrustes: rmse 0.0574152  max resid 0.2826642 
Run 30 stress 0.1736522 
Run 31 stress 0.1733041 
Run 32 stress 0.1712221 
Run 33 stress 0.1720533 
Run 34 stress 0.1733034 
Run 35 stress 0.1727863 
Run 36 stress 0.1762042 
Run 37 stress 0.1767485 
Run 38 stress 0.1710044 
Run 39 stress 0.1703045 
... Procrustes: rmse 0.05793853  max resid 0.2825322 
Run 40 stress 0.1701603 
... Procrustes: rmse 0.00243186  max resid 0.007204675 
... Similar to previous best
*** Solution reached

T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T2_rare_plot<-T2_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none")  + xlab("NMDS 1") + ylab("NMDS 2")
T2_rare_plot

Warning: Removed 4 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (geom_point).

ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1685562 
Run 1 stress 0.1748996 
Run 2 stress 0.1726844 
Run 3 stress 0.1728187 
Run 4 stress 0.1708916 
Run 5 stress 0.1723652 
Run 6 stress 0.1708747 
Run 7 stress 0.1718667 
Run 8 stress 0.1686699 
... Procrustes: rmse 0.009687857  max resid 0.05462871 
Run 9 stress 0.1684508 
... New best solution
... Procrustes: rmse 0.0119619  max resid 0.06653228 
Run 10 stress 0.1772936 
Run 11 stress 0.1687617 
... Procrustes: rmse 0.02053262  max resid 0.105457 
Run 12 stress 0.1686702 
... Procrustes: rmse 0.0151362  max resid 0.06254477 
Run 13 stress 0.1686729 
... Procrustes: rmse 0.02086917  max resid 0.1253609 
Run 14 stress 0.1792488 
Run 15 stress 0.1694632 
Run 16 stress 0.1696437 
Run 17 stress 0.170417 
Run 18 stress 0.1823142 
Run 19 stress 0.1687895 
... Procrustes: rmse 0.02300242  max resid 0.1239071 
Run 20 stress 0.1687873 
... Procrustes: rmse 0.02263201  max resid 0.1209943 
Run 21 stress 0.1723288 
Run 22 stress 0.1697892 
Run 23 stress 0.1694646 
Run 24 stress 0.1686697 
... Procrustes: rmse 0.01566615  max resid 0.06646188 
Run 25 stress 0.1720066 
Run 26 stress 0.1708448 
Run 27 stress 0.172342 
Run 28 stress 0.1708482 
Run 29 stress 0.1686102 
... Procrustes: rmse 0.01107215  max resid 0.06307271 
Run 30 stress 0.1686539 
... Procrustes: rmse 0.01829401  max resid 0.1088791 
Run 31 stress 0.1725083 
Run 32 stress 0.1729314 
Run 33 stress 0.1823118 
Run 34 stress 0.171603 
Run 35 stress 0.1684516 
... Procrustes: rmse 0.000333685  max resid 0.001749564 
... Similar to previous best
*** Solution reached

T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)

T3_rare_plot<-T3_rare$plot + theme_classic() + xlim(-0.5, 0.5) + ylim(-0.5, 0.5) + theme(legend.position = "none")  + xlab("NMDS 1") + ylab("NMDS 2")
T3_rare_plot

Warning: Removed 4 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 4 rows containing missing values (geom_point).

library(ggpubr)
ggarrange(T1_rare_plot, T2_rare_plot, T3_rare_plot, ncol = 1)

Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_ordination.pdf", width = 5, height = 10)

CAP ordination plots rarefied

t1_dist <- distance(subset_samples(ps_rare, Time=="T1"), method="bray") #get wUnifrac and save
t1_table<-as.matrix(dist(t1_dist)) #transform wUnifrac index
ord_t1_rare_cap <- capscale(t1_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T1"))))
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T1_cap.pdf")


t2_dist <- distance(subset_samples(ps_rare, Time=="T2"), method="bray") #get wUnifrac and save
t2_table<-as.matrix(dist(t2_dist)) #transform wUnifrac index
ord_t2_rare_cap <- capscale(t2_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T2"))))
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T2_cap.pdf")


#G166SG identified as outlier based on plots with it included. Removed to create plot. 
ps_rare <-  subset_samples(ps_rare, sample_names(ps_rare) != "G166SG")
t3_dist <- distance(subset_samples(ps_rare, Time=="T3"), method="bray") #get wUnifrac and save
t3_table<-as.matrix(dist(t3_dist)) #transform wUnifrac index
ord_t3_rare_cap <- capscale(t3_table ~ Herbicide, data.frame(sample_data(subset_samples(ps_rare, Time == "T3"))))
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare_cap, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_ordination_rare_T3_cap.pdf")

Pairwise adonis testing no needed becasue of insignificant gloabal test.

ps_pairwiseadonis<-function(physeq){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
pairwiseAdonis::pairwise.adonis(x = physeq_dist, factors = md_tab$Herbicide, p.adjust.m = "none", perm = 1000)
}

ps_t1<-subset_samples(ps_rare_sub, Time == "T1")
ps_t1<-prune_taxa(taxa_sums(ps_t1) > 2, ps_t1)

ps_t2<-subset_samples(ps_rare_sub, Time == "T2")
ps_t2<-prune_taxa(taxa_sums(ps_t2) > 2, ps_t2)

ps_t3<-subset_samples(ps_rare_sub, Time == "T3")
ps_t3<-prune_taxa(taxa_sums(ps_t3) > 2, ps_t3)


ps_pairwiseadonis(ps_t1)
ps_pairwiseadonis(ps_t2)
ps_pairwiseadonis(ps_t3)

Pairwise betadispr by treatment, time and mode

ps_betadispr<-function(physeq, groupingvar = "Groupingvar"){
  otu_tab<-data.frame(phyloseq::otu_table(physeq))
  md_tab<-data.frame(phyloseq::sample_data(physeq))
    if(taxa_are_rows(physeq)== T){
       physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
            else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
                mod<-vegan::betadisper(physeq_dist, md_tab[,groupingvar])
        ## Perform test
                print(anova(mod))
        ## Permutation test for F
                pmod <- vegan::permutest(mod, permutations = 1000, pairwise = TRUE)
                print(pmod)
                print(boxplot(mod))
}

permute test of dispersion

ps_betadispr(subset_samples(ps_rare_sub, Time == "T1"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T2"), groupingvar = "Mode")
ps_betadispr(subset_samples(ps_rare_sub, Time == "T3"), groupingvar = "Mode")

ps_betadispr(subset_samples(ps_rare_sub, Mode == "Chemical"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Non-Treated"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Mode == "Hand"), groupingvar = "Time")


ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Glyphosate"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Atrazine-Mesotrione"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Dicamba"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Handweeded"), groupingvar = "Time")
ps_betadispr(subset_samples(ps_rare_sub, Herbicide == "Non-Treated"), groupingvar = "Time")

ps_betadispr(ps_rare_sub, groupingvar = "Time")

Box and whisker plots of distance within group distances

#remotes::install_github("antonioggsousa/micrUBIfuns")
library(micrUBIfuns)
T1_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
T1_beta_plot <- T1_beta$plot
T1_beta_plot <- T1_beta_plot + theme_classic()+ guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T1_beta_plot

Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).

my_legend <- get_legend(T1_beta_plot)

Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_legend.pdf")

Saving 7.29 x 4.51 in image

T1_beta_plot<-T1_beta_plot+ theme(legend.position = "none") 
T1_beta_plot

Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).

T1_beta_df<- T1_beta$data
T1_betamod<-aov(formula = beta_div_value ~ group ,data = T1_beta_df)
summary(T1_betamod)

             Df Sum Sq Mean Sq F value  Pr(>F)    
group         4 0.2967 0.07419    8.83 1.1e-06 ***
Residuals   250 2.1003 0.00840                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T1_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T1_beta_df)

$group
                                        diff         lwr          upr     p adj
Dicamba-Atrazine-Mesotrione     -0.027993939 -0.07861966  0.022631781 0.5508308
Glyphosate-Atrazine-Mesotrione  -0.036412121 -0.08703784  0.014213599 0.2806138
Handweeded-Atrazine-Mesotrione  -0.058666667 -0.11176337 -0.005569963 0.0220819
Non-Treated-Atrazine-Mesotrione  0.040006061 -0.01061966  0.090631781 0.1940292
Glyphosate-Dicamba              -0.008418182 -0.05644596  0.039609594 0.9889807
Handweeded-Dicamba              -0.030672727 -0.08129845  0.019952993 0.4576855
Non-Treated-Dicamba              0.068000000  0.01997222  0.116027775 0.0012035
Handweeded-Glyphosate           -0.022254545 -0.07288027  0.028371175 0.7468226
Non-Treated-Glyphosate           0.076418182  0.02839041  0.124445957 0.0001747
Non-Treated-Handweeded           0.098672727  0.04804701  0.149298448 0.0000019

T2_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
T2_beta_plot <- T2_beta$plot
T2_beta_plot <- T2_beta_plot+ theme_classic() + theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
T2_beta_plot

Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

T2_beta_df<- T2_beta$data
T2_betamod<-aov(formula = beta_div_value ~ group ,data = T2_beta_df)
summary(T2_betamod)

             Df Sum Sq  Mean Sq F value  Pr(>F)   
group         4  0.080 0.019995   4.367 0.00201 **
Residuals   231  1.058 0.004578                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

TukeyHSD(x = T2_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T2_beta_df)

$group
                                        diff          lwr           upr     p adj
Dicamba-Atrazine-Mesotrione     -0.029188889 -0.070787897  0.0124101197 0.3048204
Glyphosate-Atrazine-Mesotrione  -0.020375758 -0.057770485  0.0170189703 0.5647456
Handweeded-Atrazine-Mesotrione   0.022206061 -0.015188667  0.0596007884 0.4780643
Non-Treated-Atrazine-Mesotrione -0.015555556 -0.054775477  0.0236643659 0.8113359
Glyphosate-Dicamba               0.008813131 -0.031069709  0.0486959716 0.9738007
Handweeded-Dicamba               0.051394949  0.011512109  0.0912777898 0.0043255
Non-Treated-Dicamba              0.013633333 -0.027965675  0.0552323419 0.8962631
Handweeded-Glyphosate            0.042581818  0.007106064  0.0780575719 0.0097763
Non-Treated-Glyphosate           0.004820202 -0.032574526  0.0422149299 0.9966025
Non-Treated-Handweeded          -0.037761616 -0.075156344 -0.0003668883 0.0465006

T3_beta<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide") 
T3_beta$plot #+ scale_color_manual(values = c("#F8766D", "#A3A500",  "#00BF7D", "#00B0F6", "#E76BF3")) + 
T3_beta_plot <- T3_beta$plot
T3_beta_plot <- T3_beta_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + ggtitle("")
T3_beta_plot <-T3_beta_plot + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.3, 0.75)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_withingroup_beta.pdf")

Saving 7.29 x 4.51 in image
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

T3_beta_df<- T3_beta$data
T3_betamod<-aov(formula = beta_div_value ~ group ,data = T3_beta_df)
summary(T3_betamod)

             Df Sum Sq  Mean Sq F value Pr(>F)
group         4 0.0672 0.016792   1.689  0.153
Residuals   252 2.5054 0.009942               

TukeyHSD(x = T3_betamod, which = "group")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ group, data = T3_beta_df)

$group
                                        diff         lwr        upr     p adj
Dicamba-Atrazine-Mesotrione      0.009452020 -0.04731207 0.06621612 0.9909320
Glyphosate-Atrazine-Mesotrione  -0.037518182 -0.08753732 0.01250096 0.2404830
Handweeded-Atrazine-Mesotrione  -0.014736364 -0.06475550 0.03528278 0.9275528
Non-Treated-Atrazine-Mesotrione -0.023409091 -0.07637301 0.02955483 0.7430333
Glyphosate-Dicamba              -0.046970202 -0.10570358 0.01176317 0.1840580
Handweeded-Dicamba              -0.024188384 -0.08292176 0.03454499 0.7896781
Non-Treated-Dicamba             -0.032861111 -0.09412180 0.02839957 0.5804264
Handweeded-Glyphosate            0.022781818 -0.02946147 0.07502511 0.7524585
Non-Treated-Glyphosate           0.014109091 -0.04096017 0.06917835 0.9554583
Non-Treated-Handweeded          -0.008672727 -0.06374199 0.04639653 0.9926689

library(ggpubr)
ggarrange(T1_beta_plot, T2_beta_plot, T3_beta_plot, ncol = 1)

Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Warning: Removed 5 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing non-finite values (stat_boxplot).
Warning: Removed 1 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta.pdf", width = 5, height = 10)

Examination of dissimliarity across time points by treatment and then again by all chemical treatments combined.

T1_beta_df$Time<-"T1"
T2_beta_df$Time<-"T2"
T3_beta_df$Time<-"T3"


beta_div_T1_T2_T3 <- rbind(T1_beta_df, T2_beta_df, T3_beta_df)

beta_anova<-function(data, Herbicide = "Herbicide"){
  df_sub<- data %>% filter(group == Herbicide)
  mod<-aov(beta_div_value ~ Time, data = df_sub)
  print(summary(mod))
  print(TukeyHSD(mod, "Time"))
  boxplot(df_sub$beta_div_value ~ df_sub$Time)
}

beta_anova(beta_div_T1_T2_T3, Herbicide = "Non-Treated")

             Df Sum Sq Mean Sq F value   Pr(>F)    
Time          2  0.234 0.11699   13.48 4.36e-06 ***
Residuals   142  1.232 0.00868                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff         lwr         upr     p adj
T2-T1 -0.09016162 -0.13451209 -0.04581115 0.0000111
T3-T1 -0.07307273 -0.11742320 -0.02872226 0.0004285
T3-T2  0.01708889 -0.02942628  0.06360405 0.6599624

beta_anova(beta_div_T1_T2_T3, Herbicide = "Handweeded")

             Df Sum Sq  Mean Sq F value Pr(>F)  
Time          2 0.0558 0.027878   4.704 0.0104 *
Residuals   152 0.9009 0.005927                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff          lwr        upr     p adj
T2-T1  0.04627273  0.009645196 0.08290026 0.0090758
T3-T1  0.03427273 -0.002354804 0.07090026 0.0719037
T3-T2 -0.01200000 -0.046747927 0.02274793 0.6928413

beta_anova(beta_div_T1_T2_T3, Herbicide = "Dicamba")

             Df Sum Sq Mean Sq F value  Pr(>F)   
Time          2 0.0733 0.03664    5.34 0.00596 **
Residuals   124 0.8507 0.00686                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff         lwr         upr     p adj
T2-T1 -0.03579495 -0.07791783 0.006327933 0.1126414
T3-T1  0.02778838 -0.01433450 0.069911266 0.2647095
T3-T2  0.06358333  0.01727131 0.109895354 0.0041226

beta_anova(beta_div_T1_T2_T3, Herbicide = "Atrazine-Mesotrione")

             Df Sum Sq Mean Sq F value Pr(>F)
Time          2 0.0292 0.01458     1.3  0.276
Residuals   153 1.7163 0.01122               
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
              diff         lwr        upr     p adj
T2-T1 -0.034600000 -0.08744615 0.01824615 0.2708607
T3-T1 -0.009657576 -0.05811807 0.03880292 0.8847098
T3-T2  0.024942424 -0.02351807 0.07340292 0.4442906

beta_anova(beta_div_T1_T2_T3, Herbicide = "Glyphosate")

             Df Sum Sq  Mean Sq F value Pr(>F)
Time          2 0.0096 0.004779   0.804  0.449
Residuals   162 0.9632 0.005945               
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = beta_div_value ~ Time, data = df_sub)

$Time
             diff         lwr        upr     p adj
T2-T1 -0.01856364 -0.05334455 0.01621728 0.4184568
T3-T1 -0.01076364 -0.04554455 0.02401728 0.7448495
T3-T2  0.00780000 -0.02698091 0.04258091 0.8564946

#regroup all chemical treatments together and rerun betadiv calcs within group. 
sample_data(ps_rare)$Mode<-sample_data(ps_rare)$Herbicide

index <- c("Dicamba", "Glyphosate", "Atrazine-Mesotrione", "Handweeded", "Non-Treated")
values <- c("Chemical", "Chemical", "Chemical", "Handweeded", "Non-Treated")

sample_data(ps_rare)$Mode<- as.factor(values[match(sample_data(ps_rare)$Mode, index)])

#+ scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 


T1_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Mode")
T1_beta_chemical_combined_plot <- T1_beta_chemical_combined$plot 
T1_beta_chemical_combined_plot<- T1_beta_chemical_combined_plot + theme_classic() + guides(color=guide_legend("Treatment")) + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T1_beta_chemical_combined_plot

Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 354 rows containing missing values (geom_point).

my_legend <- get_legend(T1_beta_chemical_combined_plot)

Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 353 rows containing missing values (geom_point).

as_ggplot(my_legend)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_beta_combined_legend.pdf")

Saving 7.29 x 4.51 in image

T1_beta_chemical_combined_plot<-T1_beta_chemical_combined_plot+ theme(legend.position = "none")
T1_beta_chemical_combined_plot

Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 353 rows containing missing values (geom_point).

T2_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Mode")
T2_beta_chemical_combined_plot <- T2_beta_chemical_combined$plot 
T2_beta_chemical_combined_plot<- T2_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T2_beta_chemical_combined_plot

Warning: Removed 373 rows containing non-finite values (stat_boxplot).
Warning: Removed 374 rows containing missing values (geom_point).

T3_beta_chemical_combined<-beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Mode")
T3_beta_chemical_combined_plot <- T3_beta_chemical_combined$plot 
T3_beta_chemical_combined_plot<- T3_beta_chemical_combined_plot + theme_classic()+ theme(legend.position = "none") + ylab("Bray-Curtis Dissimilarity") + xlab("") + theme(axis.ticks.x = element_blank(), axis.text.x = element_blank()) + ylim (0.5, 0.75) + scale_color_manual(values = c("#FFA500", "#00B0F6", "#E76BF3")) 
T3_beta_chemical_combined_plot

Warning: Removed 348 rows containing non-finite values (stat_boxplot).
Warning: Removed 349 rows containing missing values (geom_point).

ggarrange(T1_beta_chemical_combined_plot, T2_beta_chemical_combined_plot, T3_beta_chemical_combined_plot, ncol = 1)

Warning: Removed 353 rows containing non-finite values (stat_boxplot).
Warning: Removed 354 rows containing missing values (geom_point).
Warning: Removed 373 rows containing non-finite values (stat_boxplot).
Warning: Removed 375 rows containing missing values (geom_point).
Warning: Removed 348 rows containing non-finite values (stat_boxplot).
Warning: Removed 350 rows containing missing values (geom_point).

ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_within_group_beta_chemical_combined.pdf", width = 5, height = 10)

T1_beta_df_chemical_combined <- T1_beta_chemical_combined$data
T2_beta_df_chemical_combined<- T2_beta_chemical_combined$data
T3_beta_df_chemical_combined<- T3_beta_chemical_combined$data

T1_beta_df_chemical_combined$Time<-"T1"
T2_beta_df_chemical_combined$Time<-"T2"
T3_beta_df_chemical_combined$Time<-"T3"

m1<-aov(beta_div_value ~ group, data = T1_beta_df_chemical_combined)
summary(m1)
TukeyHSD(m1, "group")
boxplot(beta_div_value ~ group, data = T1_beta_df_chemical_combined)


m2<-aov(beta_div_value ~ group, data = T2_beta_df_chemical_combined)
summary(m2)
TukeyHSD(m2, "group")
boxplot(beta_div_value ~ group, data = T2_beta_df_chemical_combined)

m3<-aov(beta_div_value ~ group, data = T3_beta_df_chemical_combined)
summary(m3)
TukeyHSD(m3, "group")
boxplot(beta_div_value ~ group, data = T3_beta_df_chemical_combined)


beta_div__chemical_combined_T1_T2_T3 <- rbind(T1_beta_df_chemical_combined, T2_beta_df_chemical_combined, T3_beta_df_chemical_combined)

beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Chemical")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Hand")
beta_anova(beta_div__chemical_combined_T1_T2_T3, Herbicide = "Non-Treated")

treatment to control

plotDistances = function(p, m, s, d) {

  # calc distances
  wu = phyloseq::distance(p, m)
  wu.m = melt(as.matrix(wu))
  
  # remove self-comparisons
  wu.m = wu.m %>%
    filter(as.character(Var1) != as.character(Var2)) %>%
    mutate_if(is.factor,as.character)
  
  # get sample data (S4 error OK and expected)
  sd = data.frame(sample_data(p)) %>%
    select(s, d) %>%
    mutate_if(is.factor,as.character)
  sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
  
  # combined distances with sample data
  colnames(sd) = c("Var1", "Type1")
  wu.sd = left_join(wu.m, sd, by = "Var1")
  
  colnames(sd) = c("Var2", "Type2")
  wu.sd = left_join(wu.sd, sd, by = "Var2")
  
  #remove this line to plot all comparisons. 
  wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
  
  # plot
  ggplot(wu.sd, aes(x = Type2, y = value)) +
    theme_bw() +
    geom_point() +
    geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
    scale_color_identity() +
    facet_wrap(~ Type1, scales = "free_x") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
    ggtitle(paste0("Distance Metric = ", m))
  
}

a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T1_rare_allgroup_beta.pdf")
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T2_rare_allgroup_beta.pdf")
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_T3_rare_allgroup_beta.pdf")

library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_combined_rare_allgroup_beta.pdf", width = 7, height = 12)

Taxon abundance bar plot

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",

        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",

        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

phylumGlommed <- tax_glom(ps_rare, "Phylum")

#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_Taxon_barplot_t1.pdf")

#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t2.pdf")

#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5), 
       labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))

plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum")  + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) + 
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/ITS_Taxon_barplot_t3.pdf")

Combined herbicide and time bar plot

sample_data(ps_rare_sub)$herb_time<-paste(sample_data(ps_rare_sub)$Herbicide, sample_data(ps_rare_sub)$Time, sep = "_")
ps_rare_for_barplot <- prune_taxa(taxa_sums(ps_rare_sub) > 50, ps_rare_sub)
plot_bar(ps_rare_for_barplot, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/ITS_BarPlot_Herbicide_Time.pdf", width = 20, height = 11)

Linear modeling of abundant taxa

Tax_glom_Subset <- function (physeq, y = "taxLevel", nreturns = "Number of returns"){
   ps_1<- tax_glom(ps_rare_sub, taxrank = y )
    myTaxa <- names(sort(taxa_sums(ps_1), decreasing = TRUE)[1:nreturns])
       ps_1_sub <- prune_taxa(myTaxa, ps_1)
  return(ps_1_sub)
}

ps_rare_family_top25<-Tax_glom_Subset(physeq = ps_rare, nreturns = 25, y = "Family")
ps_rare_order_top10<-Tax_glom_Subset(physeq = ps_rare, nreturns = 10, y = "Order")

#explore top 25 taxa with plot bar
plot_bar(ps_rare_family_top25, x= "herb_time", fill = "Family") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")
plot_bar(ps_rare_order_top10, x= "herb_time", fill = "Order") + scale_fill_manual(values = col_vector) + geom_bar(stat="identity")


#write function to wrangle data prior to anova

abund_aov_wrangle <- function (physeq, y = "Tax_Level"){
  tax<-tax_table(physeq)[,y]
   meta<-data.frame(sample_data(physeq))
  counts<-data.frame(otu_table(physeq))
  rownames(counts) <- tax[,1]
  counts<-data.frame(t(counts))
   counts$Time <- meta$Time 
   counts$Herbicide <- meta$Herbicide 
  counts$Herb_time <- meta$herb_time 
  return(counts)
}            

test<-abund_aov_wrangle(ps_rare_family_top25, y = "Family")

mod_abund<-function(count_tab, IV = "Groups to be tested") {
   for(j in 1:length(unique(count_tab[,"Herbicide"]))){
         data <- count_tab %>% filter(Herbicide == unique(count_tab$Herbicide)[j])
           #change this to the number of returns from the tax_glom_subset function
   for (i in 1:25) { 
            mod <- aov(unlist(data[i]) ~ matrix(data[,IV])) 
            #sanity check
            #print(c(j, i))
   if(summary(mod)[[1]][["Pr(>F)"]][1] <= 0.05) {
    print(summary(mod))
     print(c(as.character(unique(count_tab[,"Herbicide"]))[j], names(data)[i]))
              boxplot(unlist(data[i]) ~ unlist(data[IV]), main =paste(names(data[i]), as.character(unique(count_tab[,"Herbicide"]))[j]), xlab= "Time", ylab="Abundance") 
           }
         }
      }
    }


mod_abund(test, IV = "Time")


#explore 
ps_rare_certabosidiacae <- subset_taxa(ps_rare, Family == "f:Ceratobasidiaceae")
plot_bar(ps_rare_certabosidiacae, x= "herb_time") + scale_fill_manual(values = col_vector)